home *** CD-ROM | disk | FTP | other *** search
/ C/C++ Users Group Library 1996 July / C-C++ Users Group Library July 1996.iso / vol_200 / 247_01 / bngcd.c < prev    next >
Encoding:
C/C++ Source or Header  |  1989-04-17  |  2.8 KB  |  112 lines
  1. /*
  2.  *   MIRACL Greatest Common Divisor module.
  3.  *   bngcd.c
  4.  */
  5.  
  6. #include <stdio.h>
  7. #include "miracl.h"
  8. #define abs(x)  ((x)<0? (-(x)) : (x))
  9.  
  10. /* Access global variables */
  11.  
  12. extern int   depth;    /* error tracing ..*/
  13. extern int   trace[];  /* .. mechanism    */
  14. extern bool  check;    /* overflow check  */
  15. extern small base;     /* number base     */
  16.  
  17. extern big w1;       /* workspace variables */
  18. extern big w2;
  19. extern big w0;
  20.  
  21. int gcd(x,y,z)
  22. big x;
  23. big y;
  24. big z;
  25. { /* greatest common divisor z=gcd(x,y) by Euclids  *
  26.    * method using Lehmers algorithm for big numbers */
  27.     small q,r,a,b,c,d,m;
  28.     long u,v,lq,lr;
  29.     int n;
  30.     big t;
  31.     if (ERNUM) return 0;
  32.     depth++;
  33.     trace[depth]=12;
  34.     if (TRACER) track();
  35.     copy(x,w1);
  36.     copy(y,w2);
  37.     insign(PLUS,w1);
  38.     insign(PLUS,w2);
  39.     a=b=c=d=0;
  40.     while (size(w2)!=0)
  41.     {
  42.         if (b==0)
  43.         { /* update w1 and w2 */
  44.             divide(w1,w2,w2);
  45.             t=w1,w1=w2,w2=t;   /* swap(w1,w2) */
  46.         }
  47.         else
  48.         {
  49.             premult(w1,c,z);
  50.             premult(w1,a,w1);
  51.             premult(w2,b,w0);
  52.             premult(w2,d,w2);
  53.             add(w1,w0,w1);
  54.             add(w2,z,w2);
  55.         }
  56.         if (ERNUM || size(w2)==0) break;
  57.         n=w1[0];
  58.         if (w2[0]==1)
  59.         { /* special case if w2 is now small */
  60.             r=subdiv(w1,w2[1],w1);
  61.             convert(igcd(w2[1],r),w1);
  62.             break;
  63.         }
  64.         a=1;
  65.         b=0;
  66.         c=0;
  67.         d=1;
  68.         m=w1[n]+1;
  69.         if (sizeof(long)==2*sizeof(small))
  70.         { /* use longs if they are double length */
  71.             if (n>2)
  72.             { /* squeeze out as much significance as possible */
  73.                 u=muldiv(w1[n],base,w1[n-1],m,&r);
  74.                 u=u*base+muldiv(r,base,w1[n-2],m,&r);
  75.                 v=muldiv(w2[n],base,w2[n-1],m,&r);
  76.                 v=v*base+muldiv(r,base,w2[n-2],m,&r);
  77.             }
  78.             else
  79.             {
  80.                 u=(long)base*w1[n]+w1[n-1];
  81.                 v=(long)base*w2[n]+w2[n-1];
  82.             }
  83.         }
  84.         else
  85.         {
  86.             u=muldiv(w1[n],base,w1[n-1],m,&r);
  87.             v=muldiv(w2[n],base,w2[n-1],m,&r);
  88.         }
  89.         forever
  90.         { /* work only with most significant piece */
  91.             if (((v+c)==0) || ((v+d)==0)) break;
  92.             lq=(u+a)/(v+c);
  93.             if (lq!=(u+b)/(v+d)) break;
  94.             if (lq>=MAXBASE/abs(d)) break;
  95.             q=lq;
  96.             r=a-q*c;
  97.             a=c;
  98.             c=r;
  99.             r=b-q*d;
  100.             b=d;
  101.             d=r;
  102.             lr=u-lq*v;
  103.             u=v;
  104.             v=lr;
  105.         }
  106.     }
  107.     copy(w1,z);
  108.     depth--;
  109.     return (size(w1));
  110. }
  111.  
  112.